---
title: 'Replication: html version '
output:
  html_document:
    df_print: paged
  pdf_document:
    includes:
      in_header: header.tex
    latex_engine: lualatex
    highlight: tango
    keep_tex: yes
papersize: a4
linestretch: 1
fontsize: 11pt
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


```{r}
library(tidyverse)
library(patchwork)
library(lme4)
library(sjPlot)
library(stargazer)
library(xtable)
library(extrafont)

#install.packages("extrafont")
#library(sandwich)
```


# Study 1

## Import Data

```{r}

D4 <- read.csv("data/data_study1.csv")

# More Than (MT) -> Lower Bound (LB)
# Less Than (LT) -> Upper Bound (UB)
D4$condition <- ifelse(D4$condition == "MT", "LB", "UB")

# exclude timeout decsion
D4 %>% filter(guess > 0) -> D5

```



## Table 1

```{r}

D5 %>% group_by(condition, feedback) %>% 
  summarize(
    n = n(), 
    mu = mean(mu), 
    signal = mean(signal), 
    guess = mean(guess), 
    up_bias = mean(up_bias),
    down_bias = mean(down_bias)
  ) %>% knitr::kable(digits = 3)

```


## Figure 1

```{r eval = TRUE}
D5 %>% filter(condition == "LB") %>% 
  ggplot(aes( x=signal, y = guess)) +
  geom_point(alpha = 0.2) + 
  stat_function(fun = function(x) (100+x)/2, color ="blue") + 
  stat_function(fun = function(x) (100-x)/log(100/x), color ="red") + 
  ggtitle("Lower Bound Signal") + 
  theme(plot.title=element_text(family="Times New Roman"),
        axis.title.x=element_text(family="Times New Roman"),
        axis.title.y=element_text(family="Times New Roman"), 
        legend.text=element_text(family="Times New Roman")
        ) -> p1

D5 %>% filter(condition == "UB") %>% 
  ggplot(aes( x=signal, y = guess)) +
  geom_point(alpha = 0.2) + 
  stat_function(fun = function(x) (x)/2, color ="blue") + 
  stat_function(fun = function(x) {100 - x/log(100/(100-x))}, color ="red") + 
  ggtitle("Upper Bound Signal") + 
  theme(plot.title=element_text(family="Times New Roman"),
        axis.title.x=element_text(family="Times New Roman"),
        axis.title.y=element_text(family="Times New Roman"), 
        legend.text=element_text(family="Times New Roman")
        ) -> p2

p3 <- (p1 + p2)
# 
# ggsave(filename = "pic2/fig01.pdf",
#         plot=p3,
#         width=8,
#         height=4,
#         device=cairo_pdf,
#         units = "in")

```


```{r, fig.align='center', fig.cap="Scatter plot"}
plot(p1 + p2)
#knitr::include_graphics("pic2/fig01.pdf")
```


## Figure A1

```{r eval = TRUE}

D5 %>% filter(condition == "LB", feedbackT == 0) %>% 
  ggplot(aes( x=signal, y = guess)) +
  geom_point(alpha = 0.2) + 
  stat_function(fun = function(x) (100+x)/2, color ="blue") + 
  stat_function(fun = function(x) (100-x)/log(100/x), color ="red") + 
  ggtitle("Lower Bound Signal without FB") +
  theme(plot.title=element_text(family="Times New Roman"),
        axis.title.x=element_text(family="Times New Roman"),
        axis.title.y=element_text(family="Times New Roman"), 
        legend.text=element_text(family="Times New Roman")
        ) -> p1

D5 %>% filter(condition == "LB", feedbackT == 1) %>% 
  ggplot(aes( x=signal, y = guess)) +
  geom_point(alpha = 0.2) + 
  stat_function(fun = function(x) (100+x)/2, color ="blue") + 
  stat_function(fun = function(x) (100-x)/log(100/x), color ="red") + 
  ggtitle("Lower Bound Signal with FB") + 
  theme(plot.title=element_text(family="Times New Roman"),
        axis.title.x=element_text(family="Times New Roman"),
        axis.title.y=element_text(family="Times New Roman"), 
        legend.text=element_text(family="Times New Roman")
        ) -> p2

D5 %>% filter(condition == "UB", feedbackT == 0) %>% 
  ggplot(aes( x=signal, y = guess)) +
  geom_point(alpha = 0.2) + 
  stat_function(fun = function(x) (x)/2, color ="blue") + 
  stat_function(fun = function(x) {100 - x/log(100/(100-x))}, color ="red") + 
  ggtitle("Upper Bound Signal without FB") + 
  theme(plot.title=element_text(family="Times New Roman"),
        axis.title.x=element_text(family="Times New Roman"),
        axis.title.y=element_text(family="Times New Roman"), 
        legend.text=element_text(family="Times New Roman")
        ) -> p3

D5 %>% filter(condition == "UB", feedbackT ==1) %>% 
  ggplot(aes( x=signal, y = guess)) +
  geom_point(alpha = 0.2) + 
  stat_function(fun = function(x) (x)/2, color ="blue") + 
  stat_function(fun = function(x) {100 - x/log(100/(100-x))}, color ="red") + 
  ggtitle("Upper Bound Signal with FB") + 
  theme(plot.title=element_text(family="Times New Roman"),
        axis.title.x=element_text(family="Times New Roman"),
        axis.title.y=element_text(family="Times New Roman"), 
        legend.text=element_text(family="Times New Roman")
        ) -> p4


p5 <- (p1 + p2) / (p3 + p4)

# ggsave(filename = "pic2/fig02.pdf",
#        plot=p5,
#        width=8,
#        height=6,
#        device=cairo_pdf,
#        units = "in")

```


```{r, fig.align='center', fig.cap="Scatter plot 2"}

plot( (p1 + p2) / (p3 + p4))
# knitr::include_graphics("pic2/fig02.pdf")

```


## Table 2

```{r}

D5 %>% lmer(up_bias ~ MT_dummy + (1|i), data=.) -> res1

D5 %>% lmer(up_bias ~ MT_dummy + feedback + order + round + (1|i), data=.) -> res2

```


```{r, results='asis'}
# tab_model(res1, res2,
#           transform = NULL,
#           dv.labels = c("up bias", "up bias"),
#           robust = FALSE,
#           collapse.ci = TRUE,
#           digits = 3,
#           digits.p = 3,
#           show.r2 = FALSE,
#           show.icc = FALSE,
#           show.re.var = FALSE,
#           show.ngroups = FALSE,
#           show.fstat = FALSE,
#           show.aic = TRUE,
#           show.aicc = TRUE,
#           show.dev = FALSE,
#           show.loglik = TRUE
#           )

stargazer(res1, res2, 
          covariate.labels=c("LB Dummy","feedback","order","round","Constant"),
          type = "html")
```


## Table A1

```{r, results='asis', eval = TRUE}

D5$female <- ifelse(D5$sex == "男性", 0, 1)
D5$q_score <- ifelse(D5$quiz_Q1 == "1/2", 1, 0) + 
  ifelse(D5$quiz_Q2 == "2/3", 1, 0) +
  ifelse(D5$quiz_Q3 == "1/3", 1, 0) +
  ifelse(D5$quiz_Q4 == "1/2", 1, 0) +
  ifelse(D5$quiz_Q5 == "7/12", 1, 0) +
  ifelse(D5$quiz_Q6 == "4/7", 1, 0)

D5 %>% lmer(up_bias ~ MT_dummy + (1|i), data=.) -> res1

D5 %>% lmer(up_bias ~ MT_dummy + feedback + order + round + female + q_score + (1|i), data=.) -> res2

res2

stargazer(res1, res2, 
          type = "html",
          covariate.labels=c("LB Dummy","feedback","order","round",
                            "female", "quiz score", "Constant"),
          title = 'Testing H1 and H2 (Appendix)')

```




## Table 3


```{r}

D5 %>% filter(MT_dummy == 1) %>% 
  lmer(up_bias ~ feedback + signal + order + (1|i), data=.) -> res5


D5 %>% filter(MT_dummy == 0) %>% 
  lmer(down_bias ~ feedback + signal + order + (1|i), data=.) -> res6


D5 %>% filter(MT_dummy == 1) %>% 
  lmer(up_bias ~ feedback * round + signal + order + (1|i), data=.) -> res7


D5 %>% filter(MT_dummy == 0) %>% 
  lmer(down_bias ~ feedback * round + signal + order +  (1|i), data=.) -> res8


```


```{r, results='asis'}
# tab_model(res5, res6, res7, res8, 
#           transform = NULL, 
#           dv.labels = c("up bias (MT)", "down bias (LT)", "up bias (MT)", "down bias (LT)"),
#           robust = FALSE,
#           collapse.ci = TRUE)

stargazer(res5, res6, res7, res8, 
          type = "html",
          covariate.labels=c("feedback","round", "signal","order",
                            "feedback * round", "Constant")
          )

```



## Table A2

```{r results='asis'}

D5 %>% filter(MT_dummy == 1) %>% 
  lmer(up_bias ~ feedback + signal + order + female + q_score + (1|i), data=.) -> res5


D5 %>% filter(MT_dummy == 0) %>% 
  lmer(down_bias ~ feedback + signal + order + female + q_score + (1|i), data=.) -> res6


D5 %>% filter(MT_dummy == 1) %>% 
  lmer(up_bias ~ feedback * round + signal + order + female + q_score + (1|i), data=.) -> res7


D5 %>% filter(MT_dummy == 0) %>% 
  lmer(down_bias ~ feedback * round + signal + order + female + q_score + (1|i), data=.) -> res8


stargazer(res5, res6, res7, res8, 
          type = "html",
          covariate.labels=c("feedback","round", "signal","order", "female", "quiz score",
                            "feedback * round", "Constant")
          )

```



## Table 4: Structural Estimation Result (two-type model)

```{r}

# no type change

func_ll2 <- function(par, dat, n_sub, periods){
  
  # par =  c(0.5, 0.5, 0.05)
  # dat = d1
  # n_sub = n_sub
  # periods = periods

  # sd of normal distribution
  sigma1 <- par[1]
  sigma2 <- par[2]
  
  # probability of type rational
  p <- par[3]
  
  # probability of error
  eta <- par[4]
  
  # probability of type change
  
  w <- 0

  rr <- dat$r_r
  ra <- dat$r_a
  ch <- dat$choice
  signal <- dat$signal
  MT_dummy <- dat$MT_dummy
  
  LL <- rep(NA, n_sub)
  
  for (i in 1:n_sub){
    
    lr <- rep(NA, periods[i])
    la <- rep(NA, periods[i])
    
    for (t in 1:periods[i]){

      it <- ifelse(i > 1, sum(periods[1:(i-1)]) + t, t)

      if (MT_dummy[it] == 1){
        
        # we set 100.5 as maximum and 0.5 as minimum because the range of choice is adjusted to ch +- 0.5. 
#        sig <- sigma * (100.5 - signal[it])
        sig1 <- sigma1
        sig2 <- sigma2
        
        deno_r <- pnorm(100.5, mean = rr[it], sd = sig1) - pnorm( signal[it], mean = rr[it], sd = sig1)
        deno_a <- pnorm(100.5, mean = ra[it], sd = sig2 ) - pnorm( signal[it], mean = rr[it], sd = sig2)
        
        lr[t] <- (1 - eta) * ( pnorm( ch[it] + 0.5, mean = rr[it], sd = sig1 ) - 
                               pnorm( ch[it] - 0.5, mean = rr[it], sd = sig1 )) / deno_r + eta / 100
        la[t] <- (1 - eta) * ( pnorm( ch[it] + 0.5, mean = ra[it], sd = sig2 ) - 
                               pnorm( ch[it] - 0.5, mean = ra[it], sd = sig2 )) / deno_a + eta / 100 
        
      } else {
#        sig <- sigma * (signal[it] - 0.5)        
        sig1 <- sigma1
        sig2 <- sigma2
        deno_r <- pnorm( signal[it], mean = rr[it], sd = sig1) - pnorm(0.5, mean = rr[it], sd = sig1) 
        deno_a <- pnorm( signal[it], mean = ra[it], sd = sig2) - pnorm(0.5, mean = ra[it], sd = sig2)         
        
        lr[t] <- (1 - eta) * ( pnorm( ch[it] + 0.5, mean = rr[it], sd = sig1) - 
                               pnorm( ch[it] - 0.5, mean = rr[it], sd = sig1)) / deno_r + eta / 100
        la[t] <- (1 - eta) * ( pnorm( ch[it] + 0.5, mean = ra[it], sd = sig2) - 
                               pnorm( ch[it] - 0.5, mean = ra[it], sd = sig2)) / deno_a + eta / 100
        
      }
  
    }
    
    p_a <- la %>% log %>% sum %>% exp    
    p_r <- lr %>% log %>% sum %>% exp
    
    LL[i] <- log( (p * p_r + (1-p) * p_a) )

  }
  
  # restriction
  penalty <- (w < 0) + (w > 1) + (eta < 0) + (eta > 1) + (p < 0) + (p > 1)
  
  return(LL %>% sum - 10000 * penalty)
  
}



```


```{r}

func_output <- function(res, name){
  
  v <- (- res$hessian) %>% solve 

  Est <- data.frame(
    name = name,
    par = res$par, 
    se = c( diag(v) ) %>% sqrt
  ) %>% mutate(
    z = par / se, 
    pvalue = pnorm(- abs(z)) * 2
  )
  
  LL <- res$value
  
  return(list(Est, LL))
  
}

```


```{r}

# parameter

start2 <- c(10, 0.1, 0.5, 0.05)
name2 <- c("sigma1", "sigma2",  "p", "eta")

```


```{r eval = TRUE}

D4 %>% filter(condition == "LB", feedback == 0) -> d


n_sub1 <- nrow(d) / 20
periods1 <- numeric(n_sub1)

for (i in 1:n_sub1){
  for (t in 1:20){
    periods1[i] <- periods1[i] + (d$guess[(i-1)*20 + t] > 0)
  }
}

d1 <- d %>% filter(guess > 0)



D4 %>% filter(condition == "LB", feedback == 1) -> d

n_sub2 <- nrow(d) / 20
periods2 <- numeric(n_sub2)

for (i in 1:n_sub2){
  for (t in 1:20){
    periods2[i] <- periods2[i] + (d$guess[(i-1)*20 + t] > 0)
  }
}

d2 <- d %>% filter(guess > 0)


D4 %>% filter(condition == "UB", feedback == 0) -> d

n_sub3 <- nrow(d) / 20
periods3 <- numeric(n_sub3)

for (i in 1:n_sub3){
  for (t in 1:20){
    periods3[i] <- periods3[i] + (d$guess[(i-1)*20 + t] > 0)
  }
}

d3 <- d %>% filter(guess > 0)


D4 %>% filter(condition == "UB", feedback == 1) -> d

n_sub4 <- nrow(d) / 20
periods4 <- numeric(n_sub4)

for (i in 1:n_sub4){
  for (t in 1:20){
    periods4[i] <- periods4[i] + (d$guess[(i-1)*20 + t] > 0)
  }
}

d4 <- d %>% filter(guess > 0)

```


```{r eval = FALSE}
# It takes time to execute the code in this block. 

set.seed(1234567)

suppressWarnings(
  mt0_res <- optim(par = start2, func_ll2,
        dat = d1, n_sub = n_sub1, periods = periods1,
        method = "Nelder-Mead", hessian = TRUE, control = list(fnscale = -1))
)



set.seed(1234567)

suppressWarnings(
  mt1_res <- optim(par = start2, func_ll2,
        dat = d2, n_sub = n_sub2, periods = periods2,
        method = "Nelder-Mead", hessian = TRUE, control = list(fnscale = -1))
)


set.seed(1234567)


suppressWarnings(
  lt0_res <- optim(par = start2, func_ll2,
        dat = d3, n_sub = n_sub3, periods = periods3,
        method = "Nelder-Mead", hessian = TRUE, control = list(fnscale = -1))
)


set.seed(1234567)

suppressWarnings(
  lt1_res <- optim(par = start2, func_ll2,
        dat = d4, n_sub = n_sub4, periods = periods4,
        method = "Nelder-Mead", hessian = TRUE, control = list(fnscale = -1))
)


list_res <- list(mt0_res, mt1_res, lt0_res, lt1_res)
save(list_res, file = "estimation_results/Re1.rds")
```

```{r}
load(file = "estimation_results/Re1.rds")

Re1 <- func_output(res = list_res[[1]], name = name2)
Re2 <- func_output(res = list_res[[2]], name = name2)
Re3 <- func_output(res = list_res[[3]], name = name2)
Re4 <- func_output(res = list_res[[4]], name = name2)

LL1 <- numeric(4)

Re1[[2]] -> LL1[1]
Re2[[2]] -> LL1[2]
Re3[[2]] -> LL1[3]
Re4[[2]] -> LL1[4]
```


\begin{table}[]
\caption{Estimation Result}
\centering
\begin{tabular}{lcccc}
\hline
 & \multicolumn{2}{c}{MT}             & \multicolumn{2}{c}{LT}             \\
 
 & No Feedback & Feedback & No Feedback & Feedback \\ \hline
 
$\sigma_r$ & 
`r Re1[[1]]$par[1] %>% format(digits = 4)` &  
`r Re2[[1]]$par[1] %>% format(digits = 4)` & 
`r Re3[[1]]$par[1] %>% format(digits = 4)` & 
`r Re4[[1]]$par[1] %>% format(digits = 4)` \\

      & 
(`r Re1[[1]]$se[1] %>% format(digits = 4)`) &  
(`r Re2[[1]]$se[1] %>% format(digits = 4)`) & 
(`r Re3[[1]]$se[1] %>% format(digits = 4)`) & 
(`r Re4[[1]]$se[1] %>% format(digits = 4)`) \\

&&&& \\

$\sigma_n$ & 
`r Re1[[1]]$par[2] %>% format(digits = 4)` &  
`r Re2[[1]]$par[2] %>% format(digits = 4)` & 
`r Re3[[1]]$par[2] %>% format(digits = 4)` & 
`r Re4[[1]]$par[2] %>% format(digits = 4)` \\


      & 
(`r Re1[[1]]$se[2] %>% format(digits = 4)`) &  
(`r Re2[[1]]$se[2] %>% format(digits = 4)`) & 
(`r Re3[[1]]$se[2] %>% format(digits = 4)`) & 
(`r Re4[[1]]$se[2] %>% format(digits = 4)`) \\

&&&& \\

$p_r$ & 
`r Re1[[1]]$par[3] %>% format(digits = 4)` &  
`r Re2[[1]]$par[3] %>% format(digits = 4)` & 
`r Re3[[1]]$par[3] %>% format(digits = 4)` & 
`r Re4[[1]]$par[3] %>% format(digits = 4)` \\


      & 
(`r Re1[[1]]$se[3] %>% format(digits = 4)`) &  
(`r Re2[[1]]$se[3] %>% format(digits = 4)`) & 
(`r Re3[[1]]$se[3] %>% format(digits = 4)`) & 
(`r Re4[[1]]$se[3] %>% format(digits = 4)`) \\

&&&& \\

$p_n = 1 - p_r$ & 
`r (1 - Re1[[1]]$par[3]) %>% format(digits = 4)` &  
`r (1 - Re2[[1]]$par[3]) %>% format(digits = 4)` & 
`r (1 - Re3[[1]]$par[3]) %>% format(digits = 4)` & 
`r (1 - Re4[[1]]$par[3]) %>% format(digits = 4)` \\


&&&& \\

$\eta$ & 
`r Re1[[1]]$par[4] %>% format(digits = 4)` &  
`r Re2[[1]]$par[4] %>% format(digits = 4)` & 
`r Re3[[1]]$par[4] %>% format(digits = 4)` & 
`r Re4[[1]]$par[4] %>% format(digits = 4)` \\

      & 
(`r Re1[[1]]$se[4] %>% format(digits = 4)`) &  
(`r Re2[[1]]$se[4] %>% format(digits = 4)`) & 
(`r Re3[[1]]$se[4] %>% format(digits = 4)`) & 
(`r Re4[[1]]$se[4] %>% format(digits = 4)`) \\ 

&&&& \\ \hline

Log Likelihood & 
`r Re1[[2]] %>% format(digits = 7)` &
`r Re2[[2]] %>% format(digits = 7)` & 
`r Re3[[2]] %>% format(digits = 7)` & 
`r Re4[[2]] %>% format(digits = 7)` 

\end{tabular}
\end{table}



<table>
    <caption>Estimation Result</caption>
    <thead>
        <tr>
            <th></th>
            <th colspan="2">LB</th>
            <th colspan="2">UB</th>
        </tr>
        <tr>
            <th></th>
            <th>No Feedback</th>
            <th>Feedback</th>
            <th>No Feedback</th>
            <th>Feedback</th>
        </tr>
    </thead>
    <tbody>
        <tr>
            <td>&#x03C3;<sub>r</sub></td>
            <td>`r Re1[[1]]$par[1] %>% format(digits = 4)`</td>
            <td>`r Re2[[1]]$par[1] %>% format(digits = 4)`</td>
            <td>`r Re3[[1]]$par[1] %>% format(digits = 4)`</td>
            <td>`r Re4[[1]]$par[1] %>% format(digits = 4)`</td>
        </tr>
        <tr>
            <td></td>
            <td>(`r Re1[[1]]$se[1] %>% format(digits = 4)`)</td>
            <td>(`r Re2[[1]]$se[1] %>% format(digits = 4)`)</td>
            <td>(`r Re3[[1]]$se[1] %>% format(digits = 4)`)</td>
            <td>(`r Re4[[1]]$se[1] %>% format(digits = 4)`)</td>
        </tr>
        <tr>
            <td></td>
            <td></td>
            <td></td>
            <td></td>
            <td></td>
        </tr>
        <tr>
            <td>&#x03C3;<sub>n</sub></td>
            <td>`r Re1[[1]]$par[2] %>% format(digits = 4)`</td>
            <td>`r Re2[[1]]$par[2] %>% format(digits = 4)`</td>
            <td>`r Re3[[1]]$par[2] %>% format(digits = 4)`</td>
            <td>`r Re4[[1]]$par[2] %>% format(digits = 4)`</td>
        </tr>
        <tr>
            <td></td>
            <td>(`r Re1[[1]]$se[2] %>% format(digits = 4)`)</td>
            <td>(`r Re2[[1]]$se[2] %>% format(digits = 4)`)</td>
            <td>(`r Re3[[1]]$se[2] %>% format(digits = 4)`)</td>
            <td>(`r Re4[[1]]$se[2] %>% format(digits = 4)`)</td>
        </tr>
        <tr>
            <td></td>
            <td></td>
            <td></td>
            <td></td>
            <td></td>
        </tr>
        <tr>
            <td>$p_r$</td>
            <td>`r Re1[[1]]$par[3] %>% format(digits = 4)`</td>
            <td>`r Re2[[1]]$par[3] %>% format(digits = 4)`</td>
            <td>`r Re3[[1]]$par[3] %>% format(digits = 4)`</td>
            <td>`r Re4[[1]]$par[3] %>% format(digits = 4)`</td>
        </tr>
        <tr>
            <td></td>
            <td>(`r Re1[[1]]$se[3] %>% format(digits = 4)`)</td>
            <td>(`r Re2[[1]]$se[3] %>% format(digits = 4)`)</td>
            <td>(`r Re3[[1]]$se[3] %>% format(digits = 4)`)</td>
            <td>(`r Re4[[1]]$se[3] %>% format(digits = 4)`)</td>
        </tr>
        <tr>
            <td></td>
            <td></td>
            <td></td>
            <td></td>
            <td></td>
        </tr>
        <tr>
            <td>$p_n = 1 - p_r$</td>
            <td>`r (1 - Re1[[1]]$par[3]) %>% format(digits = 4)`</td>
            <td>`r (1 - Re2[[1]]$par[3]) %>% format(digits = 4)`</td>
            <td>`r (1 - Re3[[1]]$par[3]) %>% format(digits = 4)`</td>
            <td>`r (1 - Re4[[1]]$par[3]) %>% format(digits = 4)`</td>
        </tr>
        <tr>
            <td></td>
            <td></td>
            <td></td>
            <td></td>
            <td></td>
        </tr>
        <tr>
            <td>&#x03B7;</td>
            <td>`r Re1[[1]]$par[4] %>% format(digits = 4)`</td>
            <td>`r Re2[[1]]$par[4] %>% format(digits = 4)`</td>
            <td>`r Re3[[1]]$par[4] %>% format(digits = 4)`</td>
            <td>`r Re4[[1]]$par[4] %>% format(digits = 4)`</td>
        </tr>
        <tr>
            <td></td>
            <td>(`r Re1[[1]]$se[4] %>% format(digits = 4)`)</td>
            <td>(`r Re2[[1]]$se[4] %>% format(digits = 4)`)</td>
            <td>(`r Re3[[1]]$se[4] %>% format(digits = 4)`)</td>
            <td>(`r Re4[[1]]$se[4] %>% format(digits = 4)`)</td>
        </tr>
        <tr>
            <td></td>
            <td></td>
            <td></td>
            <td></td>
            <td></td>
        </tr>
        <tr>
            <th>Log Likelihood</th>
            <td>`r Re1[[2]] %>% format(digits = 7)`</td>
            <td>`r Re2[[2]] %>% format(digits = 7)`</td>
            <td>`r Re3[[2]] %>% format(digits = 7)`</td>
            <td>`r Re4[[2]] %>% format(digits = 7)`</td>
        </tr>
    </tbody>
</table>


## Table 5: Structural Estimation Result (four-type model)


```{r}

# no type change

func_ll4 <- function(par, dat, n_sub, periods){
  
  # par =  c(0.5, 0.5, 0.05)
  # dat = d1
  # n_sub = n_sub
  # periods = periods

  # sd of rational
  sigma1 <- par[1]
  
  # sd of naive
  sigma2 <- par[2]
  
  # probability of type rational
  p <- par[3]

  # probability of type post  
  p2 <- par[4]
  
  # probability of naive
  q <- par[5]  
  
  # 1 - p - q = random
  
  # probability of error
  eta <- par[6]


  rr <- dat$r_r
  ra <- dat$r_a
  ch <- dat$choice
  signal <- dat$signal
  MT_dummy <- dat$MT_dummy
  post_d <- dat$post_d
  
  LL <- rep(NA, n_sub)
  
  for (i in 1:n_sub){
    
    lr <- rep(NA, periods[i])
    lr2 <- rep(NA, periods[i])
    la <- rep(NA, periods[i])
    la2 <- rep(NA, periods[i])
    
    for (t in 1:periods[i]){

      it <- ifelse(i > 1, sum(periods[1:(i-1)]) + t, t)
      
      if (MT_dummy[it] == 1){
        
        # we set 100.5 as maximum and 0.5 as minimum because the range of choice is adjusted to ch +- 0.5. 
#        sig <- sigma * (100.5 - signal[it])
        sig1 <- sigma1
        sig2 <- sigma2
        
        deno_r <- pnorm(100.5, mean = rr[it], sd = sig1) - pnorm( signal[it], mean = rr[it], sd = sig1)
        deno_a <- pnorm(100.5, mean = ra[it], sd = sig2 ) - pnorm( signal[it], mean = rr[it], sd = sig2)
        
        lr[t] <- (1 - eta) * ( pnorm( ch[it] + 0.5, mean = rr[it], sd = sig1 ) - 
                               pnorm( ch[it] - 0.5, mean = rr[it], sd = sig1 )) / deno_r + eta / 100
        la[t] <- (1 - eta) * ( pnorm( ch[it] + 0.5, mean = ra[it], sd = sig2 ) - 
                               pnorm( ch[it] - 0.5, mean = ra[it], sd = sig2 )) / deno_a + eta / 100 
  
        lr2[t] <- (1 - eta) * post_d[it] + eta / 100      
        la2[t] <- (1 - eta) / (100 - signal[it] + 1) + eta / 100 
        
        
      } else {
        sig1 <- sigma1
        sig2 <- sigma2
        deno_r <- pnorm( signal[it], mean = rr[it], sd = sig1) - pnorm(0.5, mean = rr[it], sd = sig1) 
        deno_a <- pnorm( signal[it], mean = ra[it], sd = sig2) - pnorm(0.5, mean = ra[it], sd = sig2)         
        
        lr[t] <- (1 - eta) * ( pnorm( ch[it] + 0.5, mean = rr[it], sd = sig1) - 
                               pnorm( ch[it] - 0.5, mean = rr[it], sd = sig1)) / deno_r + eta / 100
        la[t] <- (1 - eta) * ( pnorm( ch[it] + 0.5, mean = ra[it], sd = sig2) - 
                               pnorm( ch[it] - 0.5, mean = ra[it], sd = sig2)) / deno_a + eta / 100

        lr2[t] <- (1 - eta) * post_d[it] + eta / 100  
        la2[t] <- (1 - eta) / (signal[it]) + eta / 100 
        
      }
  
    }
    
    
    p_a <- la %>% log %>% sum %>% exp    
    p_r <- lr %>% log %>% sum %>% exp
    p_a2 <- la2 %>% log %>% sum %>% exp      
    p_r2 <- lr2 %>% log %>% sum %>% exp  
    
    # To avoid NA cases、we adjust q2 so that it does not become negative
    q2 <- max(0, 1 - p - p2 - q)
    
    LL[i] <- log( (p * p_r + p2 * p_r2 + q * p_a + q2 * p_a2) )

  }
  
  # restriction (we may delete this if we used constrOptim)
  penalty <- (eta < 0) + (eta > 1) + (p < 0) + (p > 1) + (p2 < 0) + (p2 > 1) + (q < 0) + (q > 1) + (p + p2 + q > 1)
  
  return(LL %>% sum - 1000 * penalty)
  
}



```


```{r}

# parameter

start4 <- c(10, 1, 0.2, 0.2, 0.3, 0.05)
name4 <- c("sigma1", "sigma2", "p rational", "p rational(post)", "p naive", "eta")

```


```{r eval = FALSE}

set.seed(1234567)

conv <- 1
start <- start4

while (conv != 0){
  suppressWarnings(
    mt0_res <- optim(par = start, func_ll4,
          dat = d1, n_sub = n_sub1, periods = periods1,
          method = "Nelder-Mead", hessian = TRUE, control = list(fnscale = -1))
  )
  
  start <- mt0_res$par
  conv <- mt0_res$convergence
    
}



# 2


conv <- 1
start <- start4

while (conv != 0){
  suppressWarnings(
    mt1_res <- optim(par = start, func_ll4,
          dat = d2, n_sub = n_sub2, periods = periods2,
          method = "Nelder-Mead", hessian = TRUE, control = list(fnscale = -1))
  )
  
  start <- mt1_res$par
  conv <- mt1_res$convergence
    
}


# 3

conv <- 1
start <- start4

while (conv != 0){
  suppressWarnings(
    lt0_res <- optim(par = start, func_ll4,
          dat = d3, n_sub = n_sub3, periods = periods3,
          method = "Nelder-Mead", hessian = TRUE, control = list(fnscale = -1))
  )
  
  start <- lt0_res$par
  conv <- lt0_res$convergence
    
}


# 4

conv <- 1
start <- start4

while (conv != 0){
  suppressWarnings(
    lt1_res <- optim(par = start, func_ll4,
          dat = d4, n_sub = n_sub4, periods = periods4,
          method = "Nelder-Mead", hessian = TRUE, control = list(fnscale = -1))
  )
  
  start <- lt1_res$par
  conv <- lt1_res$convergence
    
}




list_res <- list(mt0_res, mt1_res, lt0_res, lt1_res)
save(list_res, file = "estimation_results/Re2.rds")
```


```{r}
load(file = "estimation_results/Re2.rds")

Re1 <- func_output(res = list_res[[1]], name = name4) # LB, no feedback
Re2 <- func_output(res = list_res[[2]], name = name4) # LB, feedback
Re3 <- func_output(res = list_res[[3]], name = name4) # UB, no feedback
Re4 <- func_output(res = list_res[[4]], name = name4) # UB, no feedback

LL2 <- numeric(4)

Re1[[2]] -> LL2[1]
Re2[[2]] -> LL2[2]
Re3[[2]] -> LL2[3]
Re4[[2]] -> LL2[4]
```




\begin{table}[]
\caption{Estimation Result}
\centering
\begin{tabular}{lcccc}
\hline
 & \multicolumn{2}{c}{LB}             & \multicolumn{2}{c}{UB}             \\
 
 & No Feedback & Feedback & No Feedback & Feedback \\ \hline
 
$\sigma_{r, 1}$ & 
`r Re1[[1]]$par[1] %>% format(digits = 4)` &  
`r Re2[[1]]$par[1] %>% format(digits = 4)` & 
`r Re3[[1]]$par[1] %>% format(digits = 4)` & 
`r Re4[[1]]$par[1] %>% format(digits = 4)` \\

      & 
(`r Re1[[1]]$se[1] %>% format(digits = 4)`) &  
(`r Re2[[1]]$se[1] %>% format(digits = 4)`) & 
(`r Re3[[1]]$se[1] %>% format(digits = 4)`) & 
(`r Re4[[1]]$se[1] %>% format(digits = 4)`) \\

&&&& \\

$\sigma_{n, 1}$ & 
`r Re1[[1]]$par[2] %>% format(digits = 4)` &  
`r Re2[[1]]$par[2] %>% format(digits = 4)` & 
`r Re3[[1]]$par[2] %>% format(digits = 4)` & 
`r Re4[[1]]$par[2] %>% format(digits = 4)` \\


      & 
(`r Re1[[1]]$se[2] %>% format(digits = 4)`) &  
(`r Re2[[1]]$se[2] %>% format(digits = 4)`) & 
(`r Re3[[1]]$se[2] %>% format(digits = 4)`) & 
(`r Re4[[1]]$se[2] %>% format(digits = 4)`) \\

&&&& \\


$p_{r} = p_{r, 1} + p_{r, 2}$ & 
`r (Re1[[1]]$par[3] + Re1[[1]]$par[4]) %>% format(digits = 4)` &  
`r (Re2[[1]]$par[3] + Re2[[1]]$par[4]) %>% format(digits = 4)` & 
`r (Re3[[1]]$par[3] + Re3[[1]]$par[4]) %>% format(digits = 4)` & 
`r (Re4[[1]]$par[3] + Re4[[1]]$par[4]) %>% format(digits = 4)` \\


&&&& \\

$\quad p_{r, 1}$ & 
`r Re1[[1]]$par[3] %>% format(digits = 4)` &  
`r Re2[[1]]$par[3] %>% format(digits = 4)` & 
`r Re3[[1]]$par[3] %>% format(digits = 4)` & 
`r Re4[[1]]$par[3] %>% format(digits = 4)` \\


      & 
(`r Re1[[1]]$se[3] %>% format(digits = 4)`) &  
(`r Re2[[1]]$se[3] %>% format(digits = 4)`) & 
(`r Re3[[1]]$se[3] %>% format(digits = 4)`) & 
(`r Re4[[1]]$se[3] %>% format(digits = 4)`) \\

&&&& \\

$\quad p_{r, 2}$ & 
`r Re1[[1]]$par[4] %>% format(digits = 4)` &  
`r Re2[[1]]$par[4] %>% format(digits = 4)` & 
`r Re3[[1]]$par[4] %>% format(digits = 4)` & 
`r Re4[[1]]$par[4] %>% format(digits = 4)` \\

      & 
(`r Re1[[1]]$se[4] %>% format(digits = 4)`) &  
(`r Re2[[1]]$se[4] %>% format(digits = 4)`) & 
(`r Re3[[1]]$se[4] %>% format(digits = 4)`) & 
(`r Re4[[1]]$se[4] %>% format(digits = 4)`) \\ 


&&&& \\

$p_{n} = p_{n, 1} + p_{n, 2}$ & 
`r (1 - Re1[[1]]$par[3] - Re1[[1]]$par[4]) %>% format(digits = 4)` &  
`r (1 - Re2[[1]]$par[3] - Re2[[1]]$par[4]) %>% format(digits = 4)` &  
`r (1 - Re3[[1]]$par[3] - Re3[[1]]$par[4]) %>% format(digits = 4)` &  
`r (1 - Re4[[1]]$par[3] - Re4[[1]]$par[4]) %>% format(digits = 4)`  \\


&&&& \\

$\quad p_{n, 1}$ & 
`r Re1[[1]]$par[5] %>% format(digits = 4)` &  
`r Re2[[1]]$par[5] %>% format(digits = 4)` & 
`r Re3[[1]]$par[5] %>% format(digits = 4)` & 
`r Re4[[1]]$par[5] %>% format(digits = 4)` \\

      & 
(`r Re1[[1]]$se[5] %>% format(digits = 4)`) &  
(`r Re2[[1]]$se[5] %>% format(digits = 4)`) & 
(`r Re3[[1]]$se[5] %>% format(digits = 4)`) & 
(`r Re4[[1]]$se[5] %>% format(digits = 4)`) \\ 


&&&& \\

$\quad p_{n, 2} = 1 - p_{r, 1} - p_{r, 2} - p_{n, 1}$ & 
`r (1 - Re1[[1]]$par[3] - Re1[[1]]$par[4] - Re1[[1]]$par[5]) %>% format(digits = 4)` &  
`r (1 - Re2[[1]]$par[3] - Re2[[1]]$par[4] - Re2[[1]]$par[5]) %>% format(digits = 4)` &  
`r (1 - Re3[[1]]$par[3] - Re3[[1]]$par[4] - Re3[[1]]$par[5]) %>% format(digits = 4)` &  
`r (1 - Re4[[1]]$par[3] - Re4[[1]]$par[4] - Re4[[1]]$par[5]) %>% format(digits = 4)`  \\


&&&& \\

$\eta$ & 
`r Re1[[1]]$par[6] %>% format(digits = 4)` &  
`r Re2[[1]]$par[6] %>% format(digits = 4)` & 
`r Re3[[1]]$par[6] %>% format(digits = 4)` & 
`r Re4[[1]]$par[6] %>% format(digits = 4)` \\

      & 
(`r Re1[[1]]$se[6] %>% format(digits = 4)`) &  
(`r Re2[[1]]$se[6] %>% format(digits = 4)`) & 
(`r Re3[[1]]$se[6] %>% format(digits = 4)`) & 
(`r Re4[[1]]$se[6] %>% format(digits = 4)`) \\ 


&&&& \\ \hline

Log Likelihood & 
`r Re1[[2]] %>% format(digits = 7)` &
`r Re2[[2]] %>% format(digits = 7)` & 
`r Re3[[2]] %>% format(digits = 7)` & 
`r Re4[[2]] %>% format(digits = 7)` 

\end{tabular}
\end{table}



<table>
    <caption>Estimation Result</caption>
    <thead>
        <tr>
            <th></th>
            <th colspan="2">LB</th>
            <th colspan="2">UB</th>
        </tr>
        <tr>
            <th></th>
            <th>No Feedback</th>
            <th>Feedback</th>
            <th>No Feedback</th>
            <th>Feedback</th>
        </tr>
    </thead>
    <tbody>
        <tr>
            <td>&#x03C3;<sub>{r, 1}</sub></td>
            <td>`r Re1[[1]]$par[1] %>% format(digits = 4)`</td>
            <td>`r Re2[[1]]$par[1] %>% format(digits = 4)`</td>
            <td>`r Re3[[1]]$par[1] %>% format(digits = 4)`</td>
            <td>`r Re4[[1]]$par[1] %>% format(digits = 4)`</td>
        </tr>
        <tr>
            <td></td>
            <td>(`r Re1[[1]]$se[1] %>% format(digits = 4)`)</td>
            <td>(`r Re2[[1]]$se[1] %>% format(digits = 4)`)</td>
            <td>(`r Re3[[1]]$se[1] %>% format(digits = 4)`)</td>
            <td>(`r Re4[[1]]$se[1] %>% format(digits = 4)`)</td>
        </tr>
        <tr>
            <td></td>
            <td></td>
            <td></td>
            <td></td>
            <td></td>
        </tr>
        <tr>
            <td>&#x03C3;<sub>{n, 1}</sub></td>
            <td>`r Re1[[1]]$par[2] %>% format(digits = 4)`</td>
            <td>`r Re2[[1]]$par[2] %>% format(digits = 4)`</td>
            <td>`r Re3[[1]]$par[2] %>% format(digits = 4)`</td>
            <td>`r Re4[[1]]$par[2] %>% format(digits = 4)`</td>
        </tr>
        <tr>
            <td></td>
            <td>(`r Re1[[1]]$se[2] %>% format(digits = 4)`)</td>
            <td>(`r Re2[[1]]$se[2] %>% format(digits = 4)`)</td>
            <td>(`r Re3[[1]]$se[2] %>% format(digits = 4)`)</td>
            <td>(`r Re4[[1]]$se[2] %>% format(digits = 4)`)</td>
        </tr>
        <tr>
            <td></td>
            <td></td>
            <td></td>
            <td></td>
            <td></td>
        </tr>
        <tr>
            <td>$p_{r} = p_{r, 1} + p_{r, 2}$</td>
            <td>`r (Re1[[1]]$par[3] + Re1[[1]]$par[4]) %>% format(digits = 4)`</td>
            <td>`r (Re2[[1]]$par[3] + Re2[[1]]$par[4]) %>% format(digits = 4)`</td>
            <td>`r (Re3[[1]]$par[3] + Re3[[1]]$par[4]) %>% format(digits = 4)`</td>
            <td>`r (Re4[[1]]$par[3] + Re4[[1]]$par[4]) %>% format(digits = 4)`</td>
        </tr>
        <tr>
            <td></td>
            <td></td>
            <td></td>
            <td></td>
            <td></td>
        </tr>
        <tr>
            <td>&#x2003;$p_{r, 1}$</td>
            <td>`r Re1[[1]]$par[3] %>% format(digits = 4)`</td>
            <td>`r Re2[[1]]$par[3] %>% format(digits = 4)`</td>
            <td>`r Re3[[1]]$par[3] %>% format(digits = 4)`</td>
            <td>`r Re4[[1]]$par[3] %>% format(digits = 4)`</td>
        </tr>
        <tr>
            <td></td>
            <td>(`r Re1[[1]]$se[3] %>% format(digits = 4)`)</td>
            <td>(`r Re2[[1]]$se[3] %>% format(digits = 4)`)</td>
            <td>(`r Re3[[1]]$se[3] %>% format(digits = 4)`)</td>
            <td>(`r Re4[[1]]$se[3] %>% format(digits = 4)`)</td>
        </tr>
        <tr>
            <td></td>
            <td></td>
            <td></td>
            <td></td>
            <td></td>
        </tr>
        <tr>
            <td>&#x2003;$p_{r, 2}$</td>
            <td>`r Re1[[1]]$par[4] %>% format(digits = 4)`</td>
            <td>`r Re2[[1]]$par[4] %>% format(digits = 4)`</td>
            <td>`r Re3[[1]]$par[4] %>% format(digits = 4)`</td>
            <td>`r Re4[[1]]$par[4] %>% format(digits = 4)`</td>
        </tr>
        <tr>
            <td></td>
            <td>(`r Re1[[1]]$se[4] %>% format(digits = 4)`)</td>
            <td>(`r Re2[[1]]$se[4] %>% format(digits = 4)`)</td>
            <td>(`r Re3[[1]]$se[4] %>% format(digits = 4)`)</td>
            <td>(`r Re4[[1]]$se[4] %>% format(digits = 4)`)</td>
        </tr>
        <tr>
            <td></td>
            <td></td>
            <td></td>
            <td></td>
            <td></td>
        </tr>
        <tr>
            <td>$p_{n} = p_{n, 1} + p_{n, 2}$</td>
            <td>`r (1 - Re1[[1]]$par[3] - Re1[[1]]$par[4]) %>% format(digits = 4)`</td>
            <td>`r (1 - Re2[[1]]$par[3] - Re2[[1]]$par[4]) %>% format(digits = 4)`</td>
            <td>`r (1 - Re3[[1]]$par[3] - Re3[[1]]$par[4]) %>% format(digits = 4)`</td>
            <td>`r (1 - Re4[[1]]$par[3] - Re4[[1]]$par[4]) %>% format(digits = 4)`</td>
        </tr>
        <tr>
            <td></td>
            <td></td>
            <td></td>
            <td></td>
            <td></td>
        </tr>
        <tr>
            <td>&#x2003;$p_{n, 1}$</td>
            <td>`r Re1[[1]]$par[5] %>% format(digits = 4)`</td>
            <td>`r Re2[[1]]$par[5] %>% format(digits = 4)`</td>
            <td>`r Re3[[1]]$par[5] %>% format(digits = 4)`</td>
            <td>`r Re4[[1]]$par[5] %>% format(digits = 4)`</td>
        </tr>
        <tr>
            <td></td>
            <td>(`r Re1[[1]]$se[5] %>% format(digits = 4)`)</td>
            <td>(`r Re2[[1]]$se[5] %>% format(digits = 4)`)</td>
            <td>(`r Re3[[1]]$se[5] %>% format(digits = 4)`)</td>
            <td>(`r Re4[[1]]$se[5] %>% format(digits = 4)`)</td>
        </tr>
        <tr>
            <td></td>
            <td></td>
            <td></td>
            <td></td>
            <td></td>
        </tr>
        <tr>
            <td>&#x2003;$p_{n, 2} = 1 - p_{r, 1} - p_{r, 2} - p_{n, 1}$</td>
            <td>`r (1 - Re1[[1]]$par[3] - Re1[[1]]$par[4] - Re1[[1]]$par[5]) %>% format(digits = 4)`</td>
            <td>`r (1 - Re2[[1]]$par[3] - Re2[[1]]$par[4] - Re2[[1]]$par[5]) %>% format(digits = 4)`</td>
            <td>`r (1 - Re3[[1]]$par[3] - Re3[[1]]$par[4] - Re3[[1]]$par[5]) %>% format(digits = 4)`</td>
            <td>`r (1 - Re4[[1]]$par[3] - Re4[[1]]$par[4] - Re4[[1]]$par[5]) %>% format(digits = 4)`</td>
        </tr>
        <tr>
            <td></td>
            <td></td>
            <td></td>
            <td></td>
            <td></td>
        </tr>
        <tr>
            <td>$\eta$</td>
            <td>`r Re1[[1]]$par[6] %>% format(digits = 4)`</td>
            <td>`r Re2[[1]]$par[6] %>% format(digits = 4)`</td>
            <td>`r Re3[[1]]$par[6] %>% format(digits = 4)`</td>
            <td>`r Re4[[1]]$par[6] %>% format(digits = 4)`</td>
        </tr>
        <tr>
            <td></td>
            <td>(`r Re1[[1]]$se[6] %>% format(digits = 4)`)</td>
            <td>(`r Re2[[1]]$se[6] %>% format(digits = 4)`)</td>
            <td>(`r Re3[[1]]$se[6] %>% format(digits = 4)`)</td>
            <td>(`r Re4[[1]]$se[6] %>% format(digits = 4)`)</td>
        </tr>
        <tr>
            <td></td>
            <td></td>
            <td></td>
            <td></td>
            <td></td>
        </tr>
        <tr>
            <th>Log Likelihood</th>
            <td>`r Re1[[2]] %>% format(digits = 7)`</td>
            <td>`r Re2[[2]] %>% format(digits = 7)`</td>
            <td>`r Re3[[2]] %>% format(digits = 7)`</td>
            <td>`r Re4[[2]] %>% format(digits = 7)`</td>
        </tr>
    </tbody>
</table>



## The likelihood ratio test comparing the four- and two-type models

```{r}
# Do LL test

2 * (LL2[1] - LL1[1]) -> st1
2 * (LL2[2] - LL1[2]) -> st2
2 * (LL2[3] - LL1[3]) -> st3
2 * (LL2[4] - LL1[4]) -> st4

1 - pchisq(st1, 2) -> pv1
1 - pchisq(st2, 2) -> pv2
1 - pchisq(st3, 2) -> pv3
1 - pchisq(st4, 2) -> pv4

# p values of LL test
c(pv1, pv2, pv3, pv4)
```




## Figure 2



```{r}

func_ll4_post <- function(par, dat, n_sub, periods){
  
  sigma1 <- par[1]
  
  # sd of naive
  sigma2 <- par[2]
  
  # probability of type rational
  p <- par[3]

  # probability of type post  
  p2 <- par[4]
  
  # probability of naive
  q <- par[5]  
  
  # 1 - p - q = random
  
  # probability of error
  eta <- par[6]


  rr <- dat$r_r
  ra <- dat$r_a
  ch <- dat$choice
  signal <- dat$signal
  MT_dummy <- dat$MT_dummy
  post_d <- dat$post_d
  
  id <- rep(NA, n_sub)
  order <- rep(NA, n_sub)
  mt_dummy <- rep(NA, n_sub)
  post_p_r1 <- rep(NA, n_sub)
  post_p_r2 <- rep(NA, n_sub)
  post_p_n1 <- rep(NA, n_sub)
  post_p_n2 <- rep(NA, n_sub)
  
  up_bias <- rep(NA, n_sub)
  down_bias <- rep(NA, n_sub)

  
  for (i in 1:n_sub){
    
    lr <- rep(NA, periods[i])
    lr2 <- rep(NA, periods[i])
    la <- rep(NA, periods[i])
    la2 <- rep(NA, periods[i])
    
    n_up <- 0
    n_down <- 0
    
    for (t in 1:periods[i]){

      it <- ifelse(i > 1, sum(periods[1:(i-1)]) + t, t)
      n_up <- n_up + dat$up_bias[it]
      n_down <- n_down + dat$down_bias[it]
      
      if (MT_dummy[it] == 1){
        
        # we set 100.5 as maximum and 0.5 as minimum because the range of choice is adjusted to ch +- 0.5. 
#        sig <- sigma * (100.5 - signal[it])
        sig1 <- sigma1
        sig2 <- sigma2
        
        deno_r <- pnorm(100.5, mean = rr[it], sd = sig1) - pnorm( signal[it], mean = rr[it], sd = sig1)
        deno_a <- pnorm(100.5, mean = ra[it], sd = sig2 ) - pnorm( signal[it], mean = rr[it], sd = sig2)
        
        lr[t] <- (1 - eta) * ( pnorm( ch[it] + 0.5, mean = rr[it], sd = sig1 ) - 
                               pnorm( ch[it] - 0.5, mean = rr[it], sd = sig1 )) / deno_r + eta / 100
        la[t] <- (1 - eta) * ( pnorm( ch[it] + 0.5, mean = ra[it], sd = sig2 ) - 
                               pnorm( ch[it] - 0.5, mean = ra[it], sd = sig2 )) / deno_a + eta / 100 
  
        lr2[t] <- (1 - eta) * post_d[it] + eta / 100      
        la2[t] <- (1 - eta) / (100 - signal[it] + 1) + eta / 100 
        
        
      } else {
#        sig <- sigma * (signal[it] - 0.5)        
        sig1 <- sigma1
        sig2 <- sigma2
        deno_r <- pnorm( signal[it], mean = rr[it], sd = sig1) - pnorm(0.5, mean = rr[it], sd = sig1) 
        deno_a <- pnorm( signal[it], mean = ra[it], sd = sig2) - pnorm(0.5, mean = ra[it], sd = sig2)         
        
        lr[t] <- (1 - eta) * ( pnorm( ch[it] + 0.5, mean = rr[it], sd = sig1) - 
                               pnorm( ch[it] - 0.5, mean = rr[it], sd = sig1)) / deno_r + eta / 100
        la[t] <- (1 - eta) * ( pnorm( ch[it] + 0.5, mean = ra[it], sd = sig2) - 
                               pnorm( ch[it] - 0.5, mean = ra[it], sd = sig2)) / deno_a + eta / 100

        lr2[t] <- (1 - eta) * post_d[it] + eta / 100  
        la2[t] <- (1 - eta) / (signal[it]) + eta / 100 
        
      }
  
    }
    
    p_a <- la %>% log %>% sum %>% exp    
    p_r <- lr %>% log %>% sum %>% exp
    p_a2 <- la2 %>% log %>% sum %>% exp      
    p_r2 <- lr2 %>% log %>% sum %>% exp  
    
    
    deno <- p * p_r + p2 * p_r2 + q * p_a + (1 - p - p2 - q) * p_a2
    id[i] <- dat$id[it]
    order[i] <- dat$order[it]
    mt_dummy[i] <- dat$MT_dummy[it]
    post_p_r1[i] <- p * p_r / deno
    post_p_r2[i] <- p2 * p_r2 / deno
    post_p_n1[i] <- q * p_a / deno
    post_p_n2[i] <- (1 - p - p2 - q) * p_a2 / deno
    
    up_bias[i] <- n_up / periods[i]
    down_bias[i] <- n_down / periods[i]
        

  }
  
  res <- data.frame(
    id = id, 
    order = order,
    MT_dummy = mt_dummy,
    post_r1 = post_p_r1,
    post_r2 = post_p_r2,
    post_n1 = post_p_n1,
    post_n2 = post_p_n2, 
    up_bias = up_bias, 
    down_bias = down_bias
  )
  
  return(res)
  
}


```





```{r}
res1 <- list_res[[1]]
res2 <- list_res[[2]]
res3 <- list_res[[3]]
res4 <- list_res[[4]]

R1 <- func_ll4_post(par = res1$par, dat = d1, n_sub = n_sub1, periods = periods1)
R2 <- func_ll4_post(par = res2$par, dat = d2, n_sub = n_sub2, periods = periods2)
R3 <- func_ll4_post(par = res3$par, dat = d3, n_sub = n_sub3, periods = periods3)
R4 <- func_ll4_post(par = res4$par, dat = d4, n_sub = n_sub4, periods = periods4)
```


```{r}
R1 %>% mutate(
  prob_rational = (post_r1 + post_r2) / (post_r1 + post_r2 + post_n1 + post_n2), 
  prob_naive = (post_n1 + post_n2) / (post_r1 + post_r2 + post_n1 + post_n2)
) -> df1


R2 %>% mutate(
  prob_rational = (post_r1 + post_r2) / (post_r1 + post_r2 + post_n1 + post_n2), 
  prob_naive = (post_n1 + post_n2) / (post_r1 + post_r2 + post_n1 + post_n2)
) -> df2



R3 %>% mutate(
  prob_rational = (post_r1 + post_r2) / (post_r1 + post_r2 + post_n1 + post_n2), 
  prob_naive = (post_n1 + post_n2) / (post_r1 + post_r2 + post_n1 + post_n2)
) -> df3


R4 %>% mutate(
  prob_rational = (post_r1 + post_r2) / (post_r1 + post_r2 + post_n1 + post_n2), 
  prob_naive = (post_n1 + post_n2) / (post_r1 + post_r2 + post_n1 + post_n2)
) -> df4


```




```{r eval = TRUE}

df1 %>% ggplot(aes(x = up_bias, y = prob_rational)) + 
  geom_point(alpha = 0.3) +
  geom_smooth() + 
  labs(title = "LB:No Feedback") + 
  theme(plot.title=element_text(family="Times New Roman"),
        axis.title.x=element_text(family="Times New Roman"),
        axis.title.y=element_text(family="Times New Roman"), 
        legend.text=element_text(family="Times New Roman")
        ) -> p1

df2 %>% ggplot(aes(x = up_bias, y = prob_rational)) + 
  geom_point(alpha = 0.3) +
  geom_smooth()  + 
  labs(title = "LB:Feedback") + 
  theme(plot.title=element_text(family="Times New Roman"),
        axis.title.x=element_text(family="Times New Roman"),
        axis.title.y=element_text(family="Times New Roman"), 
        legend.text=element_text(family="Times New Roman")
        ) -> p2

df3 %>% ggplot(aes(x = down_bias, y = prob_rational)) + 
  geom_point(alpha = 0.3) +
  geom_smooth()  + 
  labs(title = "UB:No Feedback") + 
  theme(plot.title=element_text(family="Times New Roman"),
        axis.title.x=element_text(family="Times New Roman"),
        axis.title.y=element_text(family="Times New Roman"), 
        legend.text=element_text(family="Times New Roman")
        )-> p3

df4 %>% ggplot(aes(x = down_bias, y = prob_rational)) + 
  geom_point(alpha = 0.3) + 
  geom_smooth() + 
  labs(title = "UB:Feedback") + 
  theme(plot.title=element_text(family="Times New Roman"),
        axis.title.x=element_text(family="Times New Roman"),
        axis.title.y=element_text(family="Times New Roman"), 
        legend.text=element_text(family="Times New Roman")
        )-> p4

p5 <- (p1 + p2) / (p3 + p4)

# ggsave(filename = "pic2/fig03.pdf",
#        plot=p5,
#        width=8,
#        height=6,
#        device=cairo_pdf,
#        units = "in")

```



```{r, fig.align='center', fig.cap="Scatter plots of the posterior probability of being a rational updater and the frequency of bias"}

plot((p1 + p2) / (p3 + p4))

# knitr::include_graphics("pic2/fig03.pdf")
```



## Correlation between the probability of being a rational updater and ratio of biased choice


```{r}
cor1 <- cor.test(df1$prob_rational, df1$up_bias) # LB, no feedback
cor2 <- cor.test(df2$prob_rational, df2$up_bias) # LB, feedback
cor3 <- cor.test(df3$prob_rational, df3$down_bias) # UB, no feedback
cor4 <- cor.test(df4$prob_rational, df4$down_bias) # UB, feedback

cor1
cor2
cor3
cor4
```



## Figure 3: Scatter plots between prob. of being rational in LB and that in UB

```{r}
## LB vs UB

R1 %>% mutate(
  prob_rational = (post_r1 + post_r2) / (post_r1 + post_r2 + post_n1 + post_n2), 
  prob_naive = (post_n1 + post_n2) / (post_r1 + post_r2 + post_n1 + post_n2)
) %>% arrange(id) -> tmp1

R3 %>% mutate(
  prob_rational = (post_r1 + post_r2) / (post_r1 + post_r2 + post_n1 + post_n2), 
  prob_naive = (post_n1 + post_n2) / (post_r1 + post_r2 + post_n1 + post_n2)
) %>% arrange(id) -> tmp2

data.frame(
  p_rational_mt = tmp1$prob_rational, 
  p_rational_lt = tmp2$prob_rational 
) %>% ggplot(aes(x = p_rational_mt, y = p_rational_lt)) + 
  geom_point() + 
  labs(title = "LB vs UB:No Feedback", x ="p_rational_LB", y = "p_rational_UB") + 
  theme(plot.title=element_text(family="Times New Roman"),
        axis.title.x=element_text(family="Times New Roman"),
        axis.title.y=element_text(family="Times New Roman"), 
        legend.text=element_text(family="Times New Roman")
        ) -> p_mtlt_no_fb


data.frame(
  p_rational_mt = tmp1$prob_rational, 
  p_rational_lt = tmp2$prob_rational 
) -> df1

data.frame(
  p_rational_mt = tmp1$prob_rational, 
  p_rational_lt = tmp2$prob_rational 
) %>% mutate(
  p_rational_mt05 = ifelse(p_rational_mt > 0.5, 1, 0), 
  p_rational_lt05 = ifelse(p_rational_lt > 0.5, 1, 0), 
) -> tb01


R2 %>% mutate(
  prob_rational = (post_r1 + post_r2) / (post_r1 + post_r2 + post_n1 + post_n2), 
  prob_naive = (post_n1 + post_n2) / (post_r1 + post_r2 + post_n1 + post_n2)
) %>% arrange(id) -> tmp1

R4 %>% mutate(
  prob_rational = (post_r1 + post_r2) / (post_r1 + post_r2 + post_n1 + post_n2), 
  prob_naive = (post_n1 + post_n2) / (post_r1 + post_r2 + post_n1 + post_n2)
) %>% arrange(id) -> tmp2

data.frame(
  p_rational_mt = tmp1$prob_rational, 
  p_rational_lt = tmp2$prob_rational 
) %>% ggplot(aes(x = p_rational_mt, y = p_rational_lt)) + 
  geom_point() + 
  labs(title = "LB vs UB:Feedback", x ="p_rational_LB", y = "p_rational_UB") + 
  theme(plot.title=element_text(family="Times New Roman"),
        axis.title.x=element_text(family="Times New Roman"),
        axis.title.y=element_text(family="Times New Roman"), 
        legend.text=element_text(family="Times New Roman")
        ) -> p_mtlt_fb


data.frame(
  p_rational_mt = tmp1$prob_rational, 
  p_rational_lt = tmp2$prob_rational 
) %>% mutate(
  p_rational_mt05 = ifelse(p_rational_mt > 0.5, 1, 0), 
  p_rational_lt05 = ifelse(p_rational_lt > 0.5, 1, 0), 
) -> tb02

data.frame(
  p_rational_mt = tmp1$prob_rational, 
  p_rational_lt = tmp2$prob_rational 
) -> df2

```



```{r}
## order 1 vs order2


# w/o FB

R1 %>% mutate(
  prob_rational = (post_r1 + post_r2) / (post_r1 + post_r2 + post_n1 + post_n2), 
  prob_naive = (post_n1 + post_n2) / (post_r1 + post_r2 + post_n1 + post_n2)
) %>% filter(order == 1) -> tmp1


R3 %>% mutate(
  prob_rational = (post_r1 + post_r2) / (post_r1 + post_r2 + post_n1 + post_n2), 
  prob_naive = (post_n1 + post_n2) / (post_r1 + post_r2 + post_n1 + post_n2)
) %>% filter(order == 1) -> tmp2

Tmp1 <- rbind(tmp1, tmp2) %>% arrange(id)


R1 %>% mutate(
  prob_rational = (post_r1 + post_r2) / (post_r1 + post_r2 + post_n1 + post_n2), 
  prob_naive = (post_n1 + post_n2) / (post_r1 + post_r2 + post_n1 + post_n2)
) %>% filter(order == 2) -> tmp1


R3 %>% mutate(
  prob_rational = (post_r1 + post_r2) / (post_r1 + post_r2 + post_n1 + post_n2), 
  prob_naive = (post_n1 + post_n2) / (post_r1 + post_r2 + post_n1 + post_n2)
) %>% filter(order == 2) -> tmp2

Tmp2 <- rbind(tmp1, tmp2) %>% arrange(id)


data.frame(
  p_rational_first = Tmp1$prob_rational, 
  p_rational_second = Tmp2$prob_rational 
) %>% ggplot(aes(x = p_rational_first, y = p_rational_second)) + 
  geom_point() + 
  labs(title = "First vs Second:No Feedback", x ="p_rational_first", y = "p_rational_second") + 
  theme(plot.title=element_text(family="Times New Roman"),
        axis.title.x=element_text(family="Times New Roman"),
        axis.title.y=element_text(family="Times New Roman"), 
        legend.text=element_text(family="Times New Roman")
        ) -> p_12_no_fb


data.frame(
  p_rational_first = Tmp1$prob_rational, 
  p_rational_second = Tmp2$prob_rational 
) %>% mutate(
  p_rational_first05 = ifelse(p_rational_first > 0.5, 1, 0), 
  p_rational_second05 = ifelse(p_rational_second > 0.5, 1, 0), 
) -> tb03

data.frame(
  p_rational_first = Tmp1$prob_rational, 
  p_rational_second = Tmp2$prob_rational 
) -> df3


## w/ FB

R2 %>% mutate(
  prob_rational = (post_r1 + post_r2) / (post_r1 + post_r2 + post_n1 + post_n2), 
  prob_naive = (post_n1 + post_n2) / (post_r1 + post_r2 + post_n1 + post_n2)
) %>% filter(order == 1) -> tmp1


R4 %>% mutate(
  prob_rational = (post_r1 + post_r2) / (post_r1 + post_r2 + post_n1 + post_n2), 
  prob_naive = (post_n1 + post_n2) / (post_r1 + post_r2 + post_n1 + post_n2)
) %>% filter(order == 1) -> tmp2

Tmp1 <- rbind(tmp1, tmp2) %>% arrange(id)


R2 %>% mutate(
  prob_rational = (post_r1 + post_r2) / (post_r1 + post_r2 + post_n1 + post_n2), 
  prob_naive = (post_n1 + post_n2) / (post_r1 + post_r2 + post_n1 + post_n2)
) %>% filter(order == 2) -> tmp1


R4 %>% mutate(
  prob_rational = (post_r1 + post_r2) / (post_r1 + post_r2 + post_n1 + post_n2), 
  prob_naive = (post_n1 + post_n2) / (post_r1 + post_r2 + post_n1 + post_n2)
) %>% filter(order == 2) -> tmp2

Tmp2 <- rbind(tmp1, tmp2) %>% arrange(id)


data.frame(
  p_rational_first = Tmp1$prob_rational, 
  p_rational_second = Tmp2$prob_rational 
) %>% ggplot(aes(x = p_rational_first, y = p_rational_second)) + 
  geom_point() + 
  labs(title = "First vs Second:Feedback", x ="p_rational_first", y = "p_rational_second") + 
  theme(plot.title=element_text(family="Times New Roman"),
        axis.title.x=element_text(family="Times New Roman"),
        axis.title.y=element_text(family="Times New Roman"), 
        legend.text=element_text(family="Times New Roman")
        ) -> p_12_fb


data.frame(
  p_rational_first = Tmp1$prob_rational, 
  p_rational_second = Tmp2$prob_rational 
) %>% mutate(
  p_rational_first05 = ifelse(p_rational_first > 0.5, 1, 0), 
  p_rational_second05 = ifelse(p_rational_second > 0.5, 1, 0), 
) -> tb04

data.frame(
  p_rational_first = Tmp1$prob_rational, 
  p_rational_second = Tmp2$prob_rational 
) -> df4


p5 <- (p_mtlt_no_fb + p_mtlt_fb ) / (p_12_no_fb + p_12_fb) 
# 
# ggsave(filename = "pic2/fig04.pdf",
#        plot=p5,
#        width=8,
#        height=6,
#        device=cairo_pdf,
#        units = "in")

  
```


```{r, fig.align='center', fig.cap="Scatter plots between prob. of being rational in MT and that in LT (top two panels). Scatter plots between prob. of being rational in the first session and that in the second session (botom two panels)."}

plot( (p_mtlt_no_fb + p_mtlt_fb ) / (p_12_no_fb + p_12_fb) )

# knitr::include_graphics("pic2/fig04.pdf")

```



## Wilcoxon signed rank tests between the tasks and between the timing


```{r}

wx1 <- wilcox.test(df1$p_rational_mt, df1$p_rational_lt, paired=TRUE, exact=TRUE) 
wx2 <- wilcox.test(df2$p_rational_mt, df2$p_rational_lt, paired=TRUE, exact=FALSE) 
wx3 <- wilcox.test(df3$p_rational_first, df3$p_rational_second, paired=TRUE) 
wx4 <- wilcox.test(df4$p_rational_first, df4$p_rational_second, paired=TRUE) 

wx1 # between LB and UB, no feedback
wx2 # between LB and UB, feedback
wx3 # between order 1 and order 2, no feedback
wx4 # between order 1 and order 2, no feedback

```




# Study 2




## Import Data

```{r}

D10 <- read.csv("data/data_study2.csv")

# More Than (MT) -> Lower Bound (LB)
# Less Than (LT) -> Upper Bound (UB)
D10$condition <- ifelse(D10$condition == "MT", "LB", "UB")

# exclude timeout decsion
D10 %>% filter(guess > 0) -> D11

```


## Table 6


```{r results='asis'}

D11 %>% group_by(labo) %>% 
  summarize(
    n = n(), 
#    female = mean(female),
    mu = mean(mu), 
    signal = mean(signal), 
    guess = mean(guess), 
    up_bias = mean(up_bias),
  ) %>% knitr::kable(digits = 3)

```



## Figure A2: Scatter plot in Study 2


```{r eval = TRUE}

D11 %>% filter(condition == "LB", labo == 1) %>% 
  ggplot(aes( x=signal, y = guess)) +
  geom_point(alpha = 0.2) + 
  stat_function(fun = function(x) (100+x)/2, color ="blue") + 
  stat_function(fun = function(x) (100-x)/log(100/x), color ="red") + 
  ggtitle("Labortory") + 
  theme(plot.title=element_text(family="Times New Roman"),
        axis.title.x=element_text(family="Times New Roman"),
        axis.title.y=element_text(family="Times New Roman"), 
        legend.text=element_text(family="Times New Roman")
        ) -> p1

D11 %>% filter(condition == "LB", labo == 0) %>% 
  ggplot(aes( x=signal, y = guess)) +
  geom_point(alpha = 0.2) + 
  stat_function(fun = function(x) (100+x)/2, color ="blue") + 
  stat_function(fun = function(x) (100-x)/log(100/x), color ="red") + 
  ggtitle("Online") + 
  theme(plot.title=element_text(family="Times New Roman"),
        axis.title.x=element_text(family="Times New Roman"),
        axis.title.y=element_text(family="Times New Roman"), 
        legend.text=element_text(family="Times New Roman")
        )-> p2

p3 <- (p1 + p2)

# ggsave(filename = "pic2/fig05.pdf",
#        plot=p3,
#        width=8,
#        height=4,
#        device=cairo_pdf,
#        units = "in")

```


```{r, fig.align='center', fig.cap="Scatter plot between signal and guess (Laboratory vs Online)"}

plot(p1 + p2)
#knitr::include_graphics("pic2/fig05.pdf")

```


## Table 7


```{r}
D5 %>% filter(MT_dummy == 1, feedback == 1, order == 1) %>% select(up_bias, round, signal, female, i) -> tmp1
D11 %>% filter(labo == 0) %>% select(up_bias, round, signal, female, i) -> tmp2

tmp1$training <- 0
tmp2$training <- 1
Tmp <- rbind(tmp1, tmp2)

Tmp %>% lmer(up_bias ~ training + round + signal + (1|i), data=.) -> res01
Tmp %>% lmer(up_bias ~ training + round + signal + female +(1|i), data=.) -> res02
```



```{r}

D11 %>% lmer(up_bias ~ labo + round + signal + (1|i), data=.) -> res03
D11 %>% lmer(up_bias ~ labo + round + signal + female + (1|i), data=.) -> res04

```



```{r, results='asis', eval = TRUE}

stargazer(res01, res02, res03, res04,
          type = "html")

```


## Table A3


```{r}

D11 %>% lmer(up_bias ~ labo + round + signal + female + rule_grasp + age + econ1 + (1|i), data=.) -> res05
D11 %>% lmer(up_bias ~ labo + round + signal + female + rule_grasp + age + econ2 + (1|i), data=.) -> res06
D11 %>% lmer(up_bias ~ labo + round + signal + female + rule_grasp + age + faculty + (1|i), data=.) -> res07

```


```{r, results='asis', eval = TRUE}

stargazer(res05, res06, res07, 
          single.row=TRUE,
          type = "html")

```



## Table 8: Estimation Result (four-type model): Laboratory v.s. Online

```{r}

D10 %>% filter(labo == 0) -> d

n_sub1 <- nrow(d) / 20
periods1 <- numeric(n_sub1)

for (i in 1:n_sub1){
  for (t in 1:20){
    periods1[i] <- periods1[i] + (d$guess[(i-1)*20 + t] > 0)
  }
}

d11 <- d %>% filter(guess > 0)



D10 %>% filter(labo == 1) -> d

n_sub2 <- nrow(d) / 20
periods2 <- numeric(n_sub2)

for (i in 1:n_sub2){
  for (t in 1:20){
    periods2[i] <- periods2[i] + (d$guess[(i-1)*20 + t] > 0)
  }
}

d12 <- d %>% filter(guess > 0)


```



```{r}

# parameter

start4 <- c(10, 1, 0.2, 0.2, 0.3, 0.05)
name4 <- c("sigma1", "sigma2", "p rational", "p rational(post)", "p naive", "eta")


```


```{r eval = FALSE}

# we used the constrOptim because we encounterd the conversion problem (2022/2/11)

set.seed(1234567)

# 1 online

conv <- 1
start <- start4

while (conv != 0){
  suppressWarnings(
    online_res <- constrOptim(theta = start, func_ll4,
          ui = rbind(c(1, 0, 0, 0, 0, 0), # 1. sigma1 > 0
                 c(0, 1, 0, 0, 0, 0), # 2. sigma2 > 0
                 c(0, 0, 1, 0, 0, 0), # 3. pr1 > 0
                 c(0, 0, 1, 0, 0, 0), # 4. pr2 > 0
                 c(0, 0, 1, 0, 0, 0), # 5. pn1 > 0
                 c(0, 0, -1, -1, -1, 0), # 6 pr1 + pr2 + pn1 < 1
                 c(0, 0, 1, 0, 0, 0)), # 7. eta > 0
          ci = c(0, 0, 0, 0, 0, -1, 0),
          dat = d11, n_sub = n_sub1, periods = periods1,
          method="Nelder-Mead", mu = 1e-05, control = list(fnscale = - 1))
  )
  
  start <- online_res$par
  conv <- online_res$convergence
    
}


suppressWarnings(
  online_res <- optim(
    par = online_res$par, func_ll4, 
    dat = d11, n_sub = n_sub1, periods = periods1,
    method = "BFGS", hessian = TRUE, control = list(fnscale = -1))  
)



# 2 labo


conv <- 1
start <- start4

while (conv != 0){
  suppressWarnings(
    labo_res <- constrOptim(theta = start, func_ll4,
          ui = rbind(c(1, 0, 0, 0, 0, 0), # 1. sigma1 > 0.1 in order to avoid mis calculate the se
                 c(0, 1, 0, 0, 0, 0), # 2. sigma2 > 0.1 in order to avoid mis calculate the se
                 c(0, 0, 1, 0, 0, 0), # 3. pr1 > 0
                 c(0, 0, 1, 0, 0, 0), # 4. pr2 > 0
                 c(0, 0, 1, 0, 0, 0), # 5. pn1 > 0
                 c(0, 0, -1, -1, -1, 0), # 6 pr1 + pr2 + pn1 < 1
                 c(0, 0, 1, 0, 0, 0)), # 7. eta > 0
          ci = c(0.1, 0.1, 0, 0, 0, -1, 0),
          dat = d12, n_sub = n_sub2, periods = periods2,
          method="Nelder-Mead", mu = 1e-05, control = list(fnscale = - 1))
  )
  
  start <- labo_res$par
  conv <- labo_res$convergence
    
}


suppressWarnings(
  labo_res <- optim(
    par = labo_res$par, func_ll4, 
    dat = d12, n_sub = n_sub2, periods = periods2,
    method = "BFGS", hessian = TRUE, control = list(fnscale = -1))  
)



list_res <- list(labo_res, online_res)
save(list_res, file = "estimation_results/Re3.rds")
```


```{r}
load(file = "estimation_results/Re3.rds")

Re1 <- func_output(res = list_res[[1]], name = name4)
Re2 <- func_output(res = list_res[[2]], name = name4)

```



\begin{table}[]
\caption{Estimation Result 3: Laboratory v.s. Online}
\centering
\begin{tabular}{lcc}
\hline
 
 & Laboratory & Online  \\ \hline
 
$\sigma_{r, 1}$ & 
`r Re1[[1]]$par[1] %>% format(digits = 4)` &  
`r Re2[[1]]$par[1] %>% format(digits = 4)` \\

      & 
(`r Re1[[1]]$se[1] %>% format(digits = 4)`) &  
(`r Re2[[1]]$se[1] %>% format(digits = 4)`) \\

&& \\

$\sigma_{n, 1}$ & 
`r Re1[[1]]$par[2] %>% format(digits = 4)` &  
`r Re2[[1]]$par[2] %>% format(digits = 4)` \\


      & 
(`r Re1[[1]]$se[2] %>% format(digits = 4)`) &  
(`r Re2[[1]]$se[2] %>% format(digits = 4)`) \\

&& \\


$p_{r} = p_{r, 1} + p_{r, 2}$ & 
`r (Re1[[1]]$par[3] + Re1[[1]]$par[4]) %>% format(digits = 4)` &  
`r (Re2[[1]]$par[3] + Re2[[1]]$par[4]) %>% format(digits = 4)` \\


&& \\

$\quad p_{r, 1}$ & 
`r Re1[[1]]$par[3] %>% format(digits = 4)` &  
`r Re2[[1]]$par[3] %>% format(digits = 4)` \\


      & 
(`r Re1[[1]]$se[3] %>% format(digits = 4)`) &  
(`r Re2[[1]]$se[3] %>% format(digits = 4)`) \\

&& \\

$\quad p_{r, 2}$ & 
`r Re1[[1]]$par[4] %>% format(digits = 4)` &  
`r Re2[[1]]$par[4] %>% format(digits = 4)` \\

      & 
(`r Re1[[1]]$se[4] %>% format(digits = 4)`) &  
(`r Re2[[1]]$se[4] %>% format(digits = 4)`) \\ 


&& \\

$p_{n} = p_{n, 1} + p_{n, 2}$ & 
`r (1 - Re1[[1]]$par[3] - Re1[[1]]$par[4]) %>% format(digits = 4)` &  
`r (1 - Re2[[1]]$par[3] - Re2[[1]]$par[4]) %>% format(digits = 4)` \\


&& \\

$\quad p_{n, 1}$ & 
`r Re1[[1]]$par[5] %>% format(digits = 4)` &  
`r Re2[[1]]$par[5] %>% format(digits = 4)`  \\

      & 
(`r Re1[[1]]$se[5] %>% format(digits = 4)`) &  
(`r Re2[[1]]$se[5] %>% format(digits = 4)`)  \\ 


&& \\

$\quad p_{n, 2} = 1 - p_{r, 1} - p_{r, 2} - p_{n, 1}$ & 
`r (1 - Re1[[1]]$par[3] - Re1[[1]]$par[4] - Re1[[1]]$par[5]) %>% format(digits = 4)` &  
`r (1 - Re2[[1]]$par[3] - Re2[[1]]$par[4] - Re2[[1]]$par[5]) %>% format(digits = 4)` \\


&& \\

$\eta$ & 
`r Re1[[1]]$par[6] %>% format(digits = 4)` &  
`r Re2[[1]]$par[6] %>% format(digits = 4)` \\

      & 
(`r Re1[[1]]$se[6] %>% format(digits = 4)`) &  
(`r Re2[[1]]$se[6] %>% format(digits = 4)`)  \\ 


&& \\ \hline

Log Likelihood & 
`r Re1[[2]] %>% format(digits = 7)` &
`r Re2[[2]] %>% format(digits = 7)` 

\end{tabular}
\end{table}



<table>
    <caption>Estimation Result 3: Laboratory v.s. Online</caption>
    <thead>
        <tr>
            <th></th>
            <th>Laboratory</th>
            <th>Online</th>
        </tr>
    </thead>
    <tbody>
        <tr>
            <td>&#x03C3;<sub>{r, 1}</sub></td>
            <td>`r Re1[[1]]$par[1] %>% format(digits = 4)`</td>
            <td>`r Re2[[1]]$par[1] %>% format(digits = 4)`</td>
        </tr>
        <tr>
            <td></td>
            <td>(`r Re1[[1]]$se[1] %>% format(digits = 4)`)</td>
            <td>(`r Re2[[1]]$se[1] %>% format(digits = 4)`)</td>
        </tr>
        <tr>
            <td></td>
            <td></td>
            <td></td>
        </tr>
        <tr>
            <td>&#x03C3;<sub>{n, 1}</sub></td>
            <td>`r Re1[[1]]$par[2] %>% format(digits = 4)`</td>
            <td>`r Re2[[1]]$par[2] %>% format(digits = 4)`</td>
        </tr>
        <tr>
            <td></td>
            <td>(`r Re1[[1]]$se[2] %>% format(digits = 4)`)</td>
            <td>(`r Re2[[1]]$se[2] %>% format(digits = 4)`)</td>
        </tr>
        <tr>
            <td></td>
            <td></td>
            <td></td>
        </tr>
        <tr>
            <td>$p_{r} = p_{r, 1} + p_{r, 2}$</td>
            <td>`r (Re1[[1]]$par[3] + Re1[[1]]$par[4]) %>% format(digits = 4)`</td>
            <td>`r (Re2[[1]]$par[3] + Re2[[1]]$par[4]) %>% format(digits = 4)`</td>
        </tr>
        <tr>
            <td></td>
            <td></td>
            <td></td>
        </tr>
        <tr>
            <td>&#x2003;$p_{r, 1}$</td>
            <td>`r Re1[[1]]$par[3] %>% format(digits = 4)`</td>
            <td>`r Re2[[1]]$par[3] %>% format(digits = 4)`</td>
        </tr>
        <tr>
            <td></td>
            <td>(`r Re1[[1]]$se[3] %>% format(digits = 4)`)</td>
            <td>(`r Re2[[1]]$se[3] %>% format(digits = 4)`)</td>
        </tr>
        <tr>
            <td></td>
            <td></td>
            <td></td>
        </tr>
        <tr>
            <td>&#x2003;$p_{r, 2}$</td>
            <td>`r Re1[[1]]$par[4] %>% format(digits = 4)`</td>
            <td>`r Re2[[1]]$par[4] %>% format(digits = 4)`</td>
        </tr>
        <tr>
            <td></td>
            <td>(`r Re1[[1]]$se[4] %>% format(digits = 4)`)</td>
            <td>(`r Re2[[1]]$se[4] %>% format(digits = 4)`)</td>
        </tr>
        <tr>
            <td></td>
            <td></td>
            <td></td>
        </tr>
        <tr>
            <td>$p_{n} = p_{n, 1} + p_{n, 2}$</td>
            <td>`r (1 - Re1[[1]]$par[3] - Re1[[1]]$par[4]) %>% format(digits = 4)`</td>
            <td>`r (1 - Re2[[1]]$par[3] - Re2[[1]]$par[4]) %>% format(digits = 4)`</td>
        </tr>
        <tr>
            <td></td>
            <td></td>
            <td></td>
        </tr>
        <tr>
            <td>&#x2003;$p_{n, 1}$</td>
            <td>`r Re1[[1]]$par[5] %>% format(digits = 4)`</td>
            <td>`r Re2[[1]]$par[5] %>% format(digits = 4)`</td>
        </tr>
        <tr>
            <td></td>
            <td>(`r Re1[[1]]$se[5] %>% format(digits = 4)`)</td>
            <td>(`r Re2[[1]]$se[5] %>% format(digits = 4)`)</td>
        </tr>
        <tr>
            <td></td>
            <td></td>
            <td></td>
        </tr>
        <tr>
            <td>&#x2003;$p_{n, 2} = 1 - p_{r, 1} - p_{r, 2} - p_{n, 1}$</td>
            <td>`r (1 - Re1[[1]]$par[3] - Re1[[1]]$par[4] - Re1[[1]]$par[5]) %>% format(digits = 4)`</td>
            <td>`r (1 - Re2[[1]]$par[3] - Re2[[1]]$par[4] - Re2[[1]]$par[5]) %>% format(digits = 4)`</td>
        </tr>
        <tr>
            <td></td>
            <td></td>
            <td></td>
        </tr>
        <tr>
            <td>$\eta$</td>
            <td>`r Re1[[1]]$par[6] %>% format(digits = 4)`</td>
            <td>`r Re2[[1]]$par[6] %>% format(digits = 4)`</td>
        </tr>
        <tr>
            <td></td>
            <td>(`r Re1[[1]]$se[6] %>% format(digits = 4)`)</td>
            <td>(`r Re2[[1]]$se[6] %>% format(digits = 4)`)</td>
        </tr>
        <tr>
            <td></td>
            <td></td>
            <td></td>
        </tr>
        <tr>
            <th>Log Likelihood</th>
            <td>`r Re1[[2]] %>% format(digits = 7)`</td>
            <td>`r Re2[[2]] %>% format(digits = 7)`</td>
        </tr>
    </tbody>
</table>



